colors <- as.numeric(claslbl)+1;
}
}
return (colors);
}
PlotHeatmap<-function(mbSetObj, plotNm, smplDist, clstDist, palette, metadata,
taxrank, viewOpt, doclust, format="png", showfeatname,
appendnm, rowV=F, colV=T, var.inx=NA, border=T, width=NA, dpi=72){
mbSetObj <- .get.mbSetObj(mbSetObj);
load_pheatmap();
load_rcolorbrewer();
load_viridis();
set.seed(2805614);
#used for color pallete
variable <<- metadata;
data <- mbSetObj$dataSet$norm.phyobj;
if(mbSetObj$module.type=="mdp"){
mbSetObj$dataSet$taxa_table <- tax_table(mbSetObj$dataSet$proc.phyobj);
data <- merge_phyloseq(data, mbSetObj$dataSet$taxa_table);
}else{
taxrank <- "OTU";
}
#if more than 1500 features will be present;subset to most abundant=>1500 features.
#OTUs already in unique names;
if(ntaxa(data)>1500){
data = prune_taxa(names(sort(taxa_sums(data), TRUE))[1:1500], data);
viewOpt == "overview";
}
if(taxrank=="OTU"){
data1 <- as.matrix(otu_table(data));
rownames(data1) <- taxa_names(data);
}else{
#merging at taxonomy levels
data <- fast_tax_glom_mem(data,taxrank);
if(is.null(data)){
AddErrMsg("Errors in projecting to the selected taxanomy level!");
return(0);
}
nm <- as.character(tax_table(data)[,taxrank]);
y <- which(is.na(nm)==TRUE);
#converting NA values to unassigned
nm[y] <- "Not_Assigned";
data1 <- as.matrix(otu_table(data));
if(appendnm=="T"){
all_nm <- colnames(tax_table(data));
hg_nmindx <- which(all_nm==taxrank)-1;
if(hg_nmindx!=0){
nma <- as.character(tax_table(data)[,hg_nmindx]);
y1 <- which(is.na(nma)==TRUE);
nma[y1] <- "Not_Assigned";
nm <- paste0(nma,"_",nm);
ind <- which(nm=="Not_Assigned_Not_Assigned");
nm[ind] <- "Not_Assigned";
nm <- gsub("_Not_Assigned", "",nm, perl = TRUE);
}
}
rownames(data1) <- nm;
#all NA club together
data1 <- (t(sapply(by(data1,rownames(data1),colSums),identity)));
nm <- rownames(data1);
}
# arrange samples on the basis of slected experimental factor and using the same for annotation also
annotation <- data.frame(sample_data(data));
ind <- which(colnames(annotation)!=metadata && colnames(annotation)!="sample_id");
if(length(ind)>0){
ind1 <- ind[1];
annotation <- annotation[order(annotation[,metadata],annotation[,ind1]),];
}else{
annotation <- annotation[order(annotation[,metadata]),];
}
# remove those columns that all values are unique (continuous or non-factors)
#uniq.inx <- apply(annotation, 2, function(x){length(unique(x)) == length(x)});
#there is an additional column sample_id which need to be removed first
# get only good meta-data
good.inx <- GetDiscreteInx(annotation);
if(sum(good.inx)>0){
annotation <- annotation[,good.inx, drop=FALSE];
sam.ord <- rownames(annotation);
data1 <- data1[,sam.ord];
}else{
annotation <- NA;
}
# set up colors for heatmap
if(palette=="gbr"){
colors <- grDevices::colorRampPalette(c("green", "black", "red"), space="rgb")(256);
}else if(palette == "heat"){
colors <- grDevices::heat.colors(256);
}else if(palette == "topo"){
colors <- grDevices::topo.colors(256);
}else if(palette == "gray"){
colors <- grDevices::colorRampPalette(c("grey90", "grey10"), space="rgb")(256);
}else if(palette == "byr"){
colors <- rev(grDevices::colorRampPalette(RColorBrewer::brewer.pal(10, "RdYlBu"))(256));
}else if(palette == "viridis") {
colors <- rev(viridis::viridis(10))
}else if(palette == "plasma") {
colors <- rev(viridis::plasma(10))
}else {
colors <- rev(grDevices::colorRampPalette(RColorBrewer::brewer.pal(10, "RdBu"))(256));
}
if(showfeatname=="T"){
showfeatname<-T;
min.margin <- 360;
} else {
showfeatname<-F;
min.margin <- 200;
}
#setting the size of plot
if(is.na(width)){
minW <- 800;
myW <- ncol(data1)*20 + min.margin;
if(myW < minW){
myW <- minW;
}
w <- round(myW/72,2);
}
myH <- nrow(data1)*20 + 180;
h <- round(myH/72,2);
if(viewOpt == "overview"){
if(is.na(width)){
if(w >9.3){
w <- 9.3;
}
}
if(h > w){
h <- w;
}
}
if(border){
border.col<-"grey60";
}else{
border.col <- NA;
}
plotNm = paste(plotNm, ".", format, sep="");
mbSetObj$imgSet$heatmap<-plotNm;
if(format=="pdf"){
grDevices::pdf(file = plotNm, width=w, height=h, bg="white", onefile=FALSE);
}else{
Cairo::Cairo(file = plotNm, unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
}
# set up color schema for samples
if(palette== "gray"){
cols <- GetColorSchema(mbSetObj, T);
uniq.cols <- unique(cols);
}else{
cols <- GetColorSchema(mbSetObj, F);
uniq.cols <- unique(cols);
}
if(doclust=="T"){
rowV<-T;
}
pheatmap::pheatmap(data1,
annotation=annotation,
fontsize=8, fontsize_row=8,
clustering_distance_rows = smplDist,
clustering_distance_cols = smplDist,
clustering_method = clstDist,
show_rownames = showfeatname,
border_color = border.col,
cluster_rows = colV,
cluster_cols = rowV,
scale= "row",
color = colors
);
dev.off();
# storing for Report Generation
mbSetObj$analSet$heatmap<-data1;
mbSetObj$analSet$heatmap.dist<-smplDist;
mbSetObj$analSet$heatmap.clust<-clstDist;
mbSetObj$analSet$heat.taxalvl<-taxrank;
return(.set.mbSetObj(mbSetObj))
}
PerformMetagenomeSeqAnal<-function(mbSetObj, variable, p.lvl, shotgunid, taxrank, model){
mbSetObj <- .get.mbSetObj(mbSetObj);
load_metagenomeseq();
filt.dataphy <- mbSetObj$dataSet$filt.data;
filt.dataphy <- apply(filt.dataphy,2,as.integer);
filt.dataphy <- otu_table(filt.dataphy,taxa_are_rows =TRUE);
sample_table <- sample_data(mbSetObj$dataSet$norm.phyobj, errorIfNULL=TRUE);
filt.dataphy <- merge_phyloseq(filt.dataphy, sample_table);
taxa_names(filt.dataphy) <- rownames(mbSetObj$dataSet$filt.data);
data <- filt.dataphy;
#data<-dataSet$norm.phyobj;
cls <- as.factor(sample_data(data)[[variable]]);
lvl <- length(levels(cls));
if(mbSetObj$module.type=="mdp"){
mbSetObj$dataSet$taxa_table <- tax_table(mbSetObj$dataSet$proc.phyobj);
data <- merge_phyloseq(data, mbSetObj$dataSet$taxa_table);
}else{ #using by default names for shotgun data
taxrank <- "OTU";
}
if(taxrank!="OTU"){
#merging at taxonomy levels
data <- fast_tax_glom_mem(data, taxrank);
if(is.null(data)){
AddErrMsg("Errors in projecting to the selected taxanomy level!");
return(0);
}
nm <- as.character(tax_table(data)[,taxrank]);
#converting NA values to unassigned
nm[is.na(nm)] <- "Not_Assigned";
data1 <- as.matrix(otu_table(data));
rownames(data1) <- nm;
#all NA club together
data1 <- as.matrix(t(sapply(by(data1,rownames(data1),colSums),identity)));
data1 <- otu_table(data1, taxa_are_rows=T);
data <- merge_phyloseq(data1, sample_data(data));
nm <- taxa_names(data);
}
tree_data <<- data;
data <- phyloseq_to_metagenomeSeq(data);
data <- cumNorm(data, p=cumNormStat(data));
mod <- model.matrix(~phenoData(data)@data[,variable]);
if(model=="zigfit"){
fit <- fitZig(data, mod);
}else{
if(length(levels(cls)) > 2){
current.msg <<- paste( "More than two groups present in your experimental factor. This model can only be used with two groups.");
return(0);
}else{
fit <-fitFeatureModel(data, mod);
}
}
x <- MRfulltable(fit, number = nrow(assayData(data)$counts));
x <- x[!is.na(rownames(x)), ];
rownames(x) <- gsub(":1", "", x = rownames(x), fixed = TRUE);
x$OTUnames <- as.character(rownames(x))
if (!is.null(tax_table(data, errorIfNULL = FALSE))) {
#Attach the bacterial taxonomy to the table, if available
TAX = data.frame(tax_table(data));
TAX$OTUnames <- as.character(rownames(TAX));
res = merge(x, TAX, by = "OTUnames")
} else {
res = x;
}
# Sort and return #sighits return TRUE or FALSE
sigHits <-res$adjPvalues<=p.lvl;
de.Num <- length(which(sigHits));
if(de.Num == 0){
current.msg <<- paste( "No significant features were identified using the given p value cutoff. Please change the cutoff limit.");
}else{
current.msg <<- paste("A total of", de.Num, "significant features were identified!")
}
if(model=="ffm"){
resTable <- res[,c("pvalues","adjPvalues","logFC")];
resTable <- signif(resTable[,c(3,1,2)], digits = 5);
colnames(resTable) <- c("logFC","Pvalues","FDR");
}else{
resTable <-res[,c("pvalues","adjPvalues")];
resTable <- signif(resTable[,1:2], digits = 5);
colnames(resTable) <- c("Pvalues","FDR");
}
ord.inx <- order(resTable$Pvalues);
resTable <- resTable[ord.inx, , drop=FALSE];
write.csv(resTable, file="metageno_de_output.csv");
if(nrow(resTable) > 500){
resTable <- resTable[1:500, ];
}
mbSetObj$analSet$metagenoseq$resTable <- mbSetObj$analSet$resTable <- data.frame(resTable);
#only getting the names of DE features
diff_ft <<- rownames(resTable)[1:de.Num];
sigfeat <- rownames(resTable);
#prepare individual boxplot
box_data <- MRcounts(data);
#subset only diff. Abundant features
box_data <- box_data[sigfeat, ];
#samples in rows
box_data <- t(box_data);
box_data <- data.frame(box_data);
colnames(box_data) <- sigfeat;
claslbl <- pData(data)[ ,variable];
box_data$class <- unlist(claslbl);
mbSetObj$analSet$boxdata <- box_data;
mbSetObj$analSet$sig.count <- de.Num;
mbSetObj$analSet$anal.type <- "metagseq";
mbSetObj$analSet$var.type <- variable;
mbSetObj$analSet$metageno.taxalvl <- taxrank;
mbSetObj$analSet$id.type <- shotgunid;
return(.set.mbSetObj(mbSetObj));
}
load_metagenomeseq <- function(){
suppressMessages(library(metagenomeSeq))
}
RF.Anal <- function(mbSetObj, treeNum, tryNum, randomOn, variable, taxrank){
load_randomforest();
mbSetObj <- .get.mbSetObj(mbSetObj);
# set up random numbers
if(is.null(mbSetObj$analSet$random.seeds)){
mbSetObj$analSet$random.seeds <- GetRandomNumbers();
mbSetObj$analSet$cur.inx <- 0;
mbSetObj$analSet$rn.seed <- mbSetObj$analSet$random.seeds[1];
}
if(randomOn == -1){
rn.sd <- 123456;
}else if(randomOn == 0){ # keep current
rn.sd <- mbSetObj$analSet$rn.seed;
}else{ # random on
cur.inx <- mbSetObj$analSet$cur.inx + 1;
rn.sd <- mbSetObj$analSet$random.seeds[cur.inx];
mbSetObj$analSet$cur.inx <- cur.inx;
}
set.seed(rn.sd);
# save the seed
mbSetObj$analSet$rn.seed <- rn.sd;
if(mbSetObj$module.type=="sdp"){
taxrank<-"OTU";
data <- mbSetObj$dataSet$norm.phyobj;
data1 <- as.data.frame(t(otu_table(data)));
}else{
if(!exists("phyloseq_objs")){
phyloseq_objs <- readRDS("phyloseq_objs.RDS")
}
if(taxrank=="OTU"){
data1 <- t(phyloseq_objs$count_tables$OTU)
}else{
taxrank.inx <- which(names(phyloseq_objs$count_tables) %in% taxrank)
data1 <- t(phyloseq_objs$count_tables[[taxrank.inx]])
}
}
data.impfeat <<- data1;
cls <- sample_data(mbSetObj$dataSet$norm.phyobj)[[variable]];
variable <<- variable;
rf_out <- randomForest(data1,cls, ntree = treeNum, mtry = tryNum, importance = TRUE, proximity = TRUE);
# set up named sig table for display
impmat <- rf_out$importance;
impmat <- impmat[rev(order(impmat[,"MeanDecreaseAccuracy"])),]
sigmat <- impmat[,"MeanDecreaseAccuracy", drop=F];
sigmat <- signif(sigmat, 5);
write.csv(sigmat,file="randomforests_sigfeatures.csv");
mbSetObj$analSet$cls <- cls;
mbSetObj$analSet$rf <- rf_out;
mbSetObj$analSet$rf.sigmat <- sigmat;
return(.set.mbSetObj(mbSetObj))
}
#'@export
PlotRF.Classify<-function(mbSetObj, feature, imgName, format="png", dpi=72, width=NA){
mbSetObj <- .get.mbSetObj(mbSetObj);
imgName = paste(imgName,".", format, sep="");
mbSetObj$imgSet$rf.cls <- imgName;
if(is.na(width)){
if(feature < 5 ){
h <- feature*1.2;
w <- 9;
} else if (feature < 10){
h <- feature*1.4;
w <- 9;
} else if (feature < 15){
h <- feature/1.6;
w <- 9;
} else if (feature < 20){
h <- feature/1.8;
w <- 9;
} else if (feature < 25){
h <- feature/2;
w <- 9;
} else if (feature < 30){
h <- feature/2.2;
w <- 9;
} else if (feature < 40){
h <- feature/2.5;
w <- 9;
} else {
h <- feature/10;
w <- 9;
}
}else if(width == 0){
w <- 8;
}else{
w <- width;
}
Cairo::Cairo(file = imgName, unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
par(mar=c(4,4,3,2));
cols <- grDevices::rainbow(length(levels(mbSetObj$analSet$cls))+1);
plot(mbSetObj$analSet$rf, main="Random Forest Classification", col=cols);
legend("topright", legend = c("Overall", levels(mbSetObj$analSet$cls)), lty=2, lwd=1, col=cols);
dev.off();
return(.set.mbSetObj(mbSetObj))
}
#'@export
PlotRF.VIP<-function(mbSetObj, feature, imgName, format="png", dpi=72, width=NA){
mbSetObj <- .get.mbSetObj(mbSetObj);
imgName = paste(imgName, ".", format, sep="");
mbSetObj$imgSet$rf.imp <- imgName;
vip.score <- rev(sort(mbSetObj$analSet$rf$importance[,"MeanDecreaseAccuracy"]));
cls <- sample_data(mbSetObj$dataSet$norm.phyobj)[[variable]];
cls.length <- length(levels(cls))
cls.width <- cls.length/2
if(is.na(width)){
if(feature < 5 ){
h <- feature*1.2;
} else if (feature < 10){
h <- feature*1.4;
} else if (feature < 15){
h <- feature/1.6;
} else if (feature < 20){
h <- feature/1.8;
} else if (feature < 25){
h <- feature/2;
} else if (feature < 30){
h <- feature/2.2;
} else if (feature < 40){
h <- feature/2.5;
} else {
h <- feature/10;
}
if(cls.length < 5){
w <- 9.25;
} else if(cls.length < 10){
w <- 11.5;
} else if(cls.length < 15){
w <- 12.5;
}else{
w <- 15;
}
}else if(width == 0){
w <- 8;
}else{
w <- width;
}
Cairo::Cairo(file = imgName,  unit="in", dpi=dpi, width=w, height=h, type=format, bg="white");
PlotImpVar(mbSetObj, vip.score,"MeanDecreaseAccuracy", feature);
dev.off();
return(.set.mbSetObj(mbSetObj))
}
PlotImpVar <- function(mbSetObj, imp.vec, xlbl, feature, color.BW=FALSE){
mbSetObj <- .get.mbSetObj(mbSetObj);
cls.len <- length(levels(mbSetObj$analSet$cls));
if(cls.len == 2){
rt.mrg <- 5;
}else if(cls.len == 3){
rt.mrg <- 6;
}else if(cls.len == 4){
rt.mrg <- 7;
}else if(cls.len == 5){
rt.mrg <- 8;
}else if(cls.len == 6){
rt.mrg <- 9;
}else{
rt.mrg <- 11;
}
op <- par(mar=c(5,9,2,rt.mrg)); # set right side margin with the number of class
feat.num = feature;
if(feat.num <= 0){
feat.num = 15;
}
if(feat.num > length(imp.vec)){
feat.num <- length(imp.vec);
}
# first get the top subset
imp.vec <- rev(sort(imp.vec))[1:feat.num];
# reverser the order for display
imp.vec <- sort(imp.vec);
# as data should already be normalized, use mean/median should be the same
# mns is a list contains means of all vars at each level
# conver the list into a matrix with each row contains var averages across different lvls
data1 <- data.impfeat;
mns <- by(data1[, names(imp.vec)], mbSetObj$analSet$cls,
function(x){ # inner function note, by send a subset of dataframe
apply(x, 2, mean, trim=0.1)
});
mns <- t(matrix(unlist(mns), ncol=feat.num, byrow=TRUE));
vip.nms <- names(imp.vec);
names(imp.vec) <- NULL;
# modified for B/W color
dotcolor <- ifelse(color.BW, "darkgrey", "#585855");
dotchart(imp.vec, bg=dotcolor, xlab= xlbl, cex=1.35);
mtext(side=2, at=1:feat.num, vip.nms, las=2, line=1, cex=1.1)
axis.lims <- par("usr"); # x1, x2, y1 ,y2
# get character width
shift <- 2*par("cxy")[1];
lgd.x <- axis.lims[2] + shift;
x <- rep(lgd.x, feat.num);
y <- 1:feat.num;
par(xpd=T);
load_rcolorbrewer();
nc <- ncol(mns);
# modified for B/W color
colorpalette <- ifelse(color.BW, "Greys", "RdYlBu");
col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(10, colorpalette))(nc); # set colors for each class
if(color.BW) col <- rev(col);
# calculate background
bg <- matrix("", nrow(mns), nc);
for (m in 1:nrow(mns)){
bg[m,] <- (col[nc:1])[rank(mns[m,])];
}
cls.lbl <- levels(mbSetObj$analSet$cls);
for (n in 1:ncol(mns)){
points(x,y, bty="n", pch=22, bg=bg[,n], cex=3);
# now add label
text(x[1], axis.lims[4], cls.lbl[n], srt=45, adj=c(0.2,0.5), cex=1.1);
# shift x, note, this is good for current size
x <- x + shift/1.25;
}
# now add color key, padding with more intermediate colors for contiuous band
col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(25, colorpalette))(50)
if(color.BW) col <- rev(col);
nc <- length(col);
x <- rep(x[1] + shift, nc);
shifty <- (axis.lims[4]-axis.lims[3])/3;
starty <- axis.lims[3] + shifty;
endy <- axis.lims[3] + 2*shifty;
y <- seq(from = starty, to = endy, length = nc);
points(x,y, bty="n", pch=15, col=rev(col), cex=2);
text(x[1], endy+shifty/8, "High", cex=1.1);
text(x[1], starty-shifty/8, "Low", cex=1.1);
par(op);
}
GetRandomNumbers <- function(){
rm(.Random.seed);
runif(1);
return(123456);
}
